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^ ', Abstract We investigate within density functional theory various physical prop- 

el ■ erties of the zero-temperature unitary Fermi gas which critically depend on the 

bX^ presence of a dispersive gradient term in the equation of state. First, we consider 

the unitary Fermi superfluid gas confined to a semi-infinite domain and calculate 
^ ■ analytically its density profile and surface tension. Then we study the quadrupole 

modes of the superfluid system under harmonic confinement finding a reliable 
analytical formula for the oscillation frequency, which reduces to the familiar 
Thomas-Fermi one in the limit of a large number of atoms. Finally, we discuss 
"j^ ■ the formation and propagation of dispersive shock waves in the collision between 

two resonant fermionic clouds, and compare our findings with recent experimental 
results. 

■ PACS numbers: 05.30.Fk, 03.75.Ss, 67.85.-d 

o: 



1 Introduction 

In the last years the crossover from the weakly paired Bardeen-Cooper-Schrieffer 
(BCS) state to the Bose-Einstein condensate (BEC) of molecular dimers with 
ultra-cold two-hyperfine-components Fermi vapors atoms has been investigated 
■ by several experimental and theoretical groups ' . When the densities of the two 



spin components are equal, and when the gas is dilute so that the range of the 
inter-atomic potential is much smaller than the inter-particle distance, then the in- 
I teraction effects are described by only one parameter: the s-wave scattering length, 

whose sign determines the character of the gas. Fano-Feshbach resonances can be 
used to change the value and the sign of the scattering length, simply by tuning 



' Dipartimento di Fisica e Astronomia "Galileo Galilei" and CNISM, Universita di Padova, Via 
Marzolo 8, 35122 Padova, Italy 



■ -CNR-IOM Democritos, via Bonomea, 265 - 34136 Trieste, Italy 

E-mail; francesco.ancilotto@unipd.it, luca.salasnich@unipd.it, flavio.toigo@unipd.it 



2 



an external magnetic field. At resonance the scattering length diverges so that the 
gas displays a very peculiar character, being at the same time dilute and strongly 
interacting. In this regime all scales associated with interactions disappear from 
the problem and the energy of the system is expected to be proportional to that of 
a non interacting fermions system. This is called the unitary regime^. 



Recently it has been remarked^ that the superfluid unitary Fermi gas, char- 
acterized by a divergent s-wave scattering length', can be efficiently described at 
zero temperature by phenomenological density functional theory. Indeed, different 
theoretical groups have proposed various density functionals. For example Bulgac 
and Yu have introduced a superfluid density functional based on a Bogoliubov- 
de Gennes approach to superfluid fermions Papenbrock and Bhattacharyya— 
have instead proposed a Kohn-Sham density functional with an effective mass to 
take into account nonlocality effects. Here we adopt instead the extended Thomas- 
Fermi functional of the unitary Fermi gas that we have proposed few years ago^-. 
The total energy in the extended Thomas-Fermi functional contains a term pro- 
portional to the kinetic energy of a uniform non interacting gas of fermions with 
number density n(r), plus a gradient correction of the form A/i^/(8m)(V«/n)^, 
originally introduced by von Weizsacker to treat surface effects in nuclei^, and 
then extensively applied to study electrons^, showing good agreement with Kohn- 
Sham calculations. In the context of the BCS-BEC crossover, the presence of the 
gradient term (and its actual weight) is a debated issue ^"'■''■'^•'^■''*-'^ '^. The main 
advantage of using such a functional is that, as it depends only on a single func- 
tion of the coordinates, i.e. the order parameter, can be used to study systems with 
quite large number of particles A^. Other functionals, which are based instead on 
single-particle orbitals, require self-consistent calculations with a numerical load 
rapidly increasing with A'^. 



In the last years we have successfully applied our extended Thomas-Fermi 
density functional and its time-dependent version^ to investigate density pro- 
files^''^, collective excitations'^, Josephson effect'^ and shock waves '^•^*' of the 
unitary Fermi gas. In addition, the collective modes of our density functional have 
been used to study the low-temperature thermodynamics of the unitary Fermi 
gas (superfluid fraction, first sound and second sound) and also the viscosity- 
entropy ratio of the unitary Fermi gas from zero-temperature elementary excita- 
tions^. 



For the superfluid unitary Fermi gas one expects the coexistence of dispersive 
and dissipative terms in the equation of stat e ^"i^'i^^i^^ . However, at zero tempera- 
ture only dispersive term survive 20 24 ^j^j^ paper we investigate various physical 
quantities of the unitary Fermi gas, which depend on the presence of a dispersive 
gradient term in the zero-temperature equation of state. In the fist part we investi- 
gate static properties like density profiles and surface tension. In the second part 
we study dynamical properties like collective oscillations and dispersive shock 
waves due to the collision between two fermionic clouds. 
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2 Extended Thomas-Fermi density functional 

The energy density of a uniform Fermi gas at unitarity depends on the constant 
density n as follows' 

^unif{n) = ^\^^fin^fl'n"' (1) 

where ^ ~ 0.4 is a universal parameter of the Fermi gas^. In the presence of an 
external trapping potential U (r) the Hohenberg-Kohn theorem ensures that the 
ground-state local density «(r) of the system can be obtained by minimizing the 
energy functional 

E[n] =F[n] + j U{r) «(r) d^r (2) 

where F[n] is a (generally unknown) energy functional of the internal energy, 
which is independent on U (r). In our extended Thomas-Fermi approach we choose 

F[n\ = j S[nyn) (3) 

where the local energy density S'{n, Vn) is given by 

^(n,Vn)=^„„,y(n)+A — (4) 
8m n 

In this expression, which is the equation of state (internal energy density) of the 
inhomogeneous system, the first term is the Thomas-Fermi-like term describing 
the uniform system, while the second term is the von Weizsacker gradient correc- 
tion^. We have obtained the value A ~ 0.25 by fitting accurate Monte Carlo results 
for the energy of fermions confined in a spherical harmonic trap close to unitary 
conditions'*^. 



3 Semi-infinite domain: density profile and surface tension 

We put in evidence the effects of the gradient term of our density functional by 
considering the external potential 

jj, s f +°° for z < 
^W = | forz>0 

acting on the Fermi superfluid. This implies that the superfluid density «(r) must 

go to zero at the boundary z = 0, while it becomes a constant n = ^ ( - ) ^ far 

from the boundary. Since in the present calculations we fix the chemical potential 
p. rather than the total number of fermions A', it is useful to introduce the zero- 
temperature grand potential energy functional of the unitary Fermi gas 



Q= [^{nyn) + U {r)n{r) - jXn{r)] dh 



(6) 
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It is convenient to introduce the characteristic length of the system 



(37r2n)i/3 



(7) 



and to rewrite the density in terms of a function of the adimensional variable ^ 

zjls as 



The grand-potential ^ then becomes 



Q = A Is nil 



/'(C)' + 7/(0'°/' -/(C)' ^c 



(8) 



(9) 



where /'(C) = df /d^ . A is the area in the {x,y) plane and L is the length in the z 
direction. The function /(C) minimizing the grand-potential obeys the equation: 



-/"(C)+/(C)'/' = /(C) 



(10) 



with the boundary conditions /(O) = /'(+°°) = and f{+°°) = 1 ■ The first 
integral of the system is 



^ = ^/'(C)' + ^/(C)'-^/(C)'°/' 

By using the appropriate boundary conditions we find K = ^ and /'(O) 
is then straightforward to get the integral equation 



^= / 

Jo 



df 



(11) 



. It 



(12) 



which gives implicitly the profile function /(C)- 

In Fig.[T]we plot the scaled density profile /(C) obtained from the numerical 

integration (solid line) of Eq. ( fT2] i. together with the function /(C) = tanh (y^lQ 
(dot-dashed line), which provides an excellent overall approximation to it. 

In the limit L/l^ — > +°° the grand potential energy Q. /A is divergent due to the 
asymptotic energy 



^asy 2 /7-T"-4/3/-._2\l/3 Z"^/'' 



A 5 2m Jo 



dC 



where /(^) = 1 and f{C) = 0. The surface tension a is defined as 



(13) 



(14) 



in the limit L/Z, 



-oo. I.e. 



a = VA^n'/'(37r2)V3. 



2m 



f 

Jo 



/'(C)' + ^/(C)'°/'-/(C)' + ^j 



d^ (15) 
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Fig. 1 Scaled density profile f{Q of the unitary Fermi gas in the semi-infinite domain. 
Solid line: numerical integration of Eq. i fTlt . Dot-dashed line: analytical interpolation f{Q = 



Eq. (fTTT i with A' = 1/5, and dl^ = df/f'{Q, helps us to rewrite the formula of the 
surface tension as 



In Fig. |2] we plot the surface tension a as a function of the von Weizsacker 
gradient coefficient A setting ^ = 0.4 By using A = 1/4, which is a reasonable 

estimate^i^ consistent with Monte Carlo results^' --, one finds a = O.SOn"^/^^ 
This value is of the same order of magnitude of previous microscopic determina- 
tions of the surface tension a based on different theoretical approaches ^^'^^'^^ and 
different boundary conditions at the interface. 



4 Extended superfluid hydrodynamics 

Let us now consider dynamical properties of the unitary Fermi gas. The starting 
point are the zero-temperature equations of superfluid hydrodynamics-^, given by 





(16) 





3- + V.(«v) = 



(17) 



(18) 
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Fig. 2 Surface tension c as a function of the von Weizsiicker gradient coefficient 1 . d is in units 

oi,tl^h^/{2m). 

where where n{r,t) is the time-dependent scalar density field and v(r,f) is the 
time-dependent vector velocity field. In the case of the unitary Fermi gas the local 
energy density can be written by using Eq. (|4]l. If /I = 0, then Eqs. (117b and dlSII 
reproduce by construction the familiar equations of superfluid hydrodynamics-^. 

For a fermionic superfluid the velocity is irrotational, i.e. V A v = 0, and the 
circulation is quantized, i.e. 



where / is an integer quantum number. This means that the velocity field of the 
unitary Fermi gas can be written as 



where 0(r,f) is the phase of the Cooper pairs i-^". Eqs. ^ and OHl can be inter- 
preted as the Euler-Lagrange equations of the following action functional 




(19) 



v(r>O = :r-V0(r,O 



(20) 




which depends on the local density n(r,f) and the phase 0(r,f). In the next sec- 
tion we use this action functional (121b to investigate zero-temperature collective 
oscillations of the unitary Fermi gas trapped by an external potential. 
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Fig. 3 Quadrupole frequency of the unitary Fermi gas witfi A' atoms under harmonic con- 
finement of frequency (O. Three different values of the gradient coefficient X. For A = (TF 
limit): Q,2 = \/2tB. Filled squares are the numerical results, while the curves are obtained with 
the analytical formula, Eq. yO). 



5 Quadrupole oscillations 

In the case of spherically-symmetric harmonic confinement, i.e. 

C/(r) = ^mw^r^ , (22) 

we have numerically studied— the collective modes of the unitary Fermi gas for 
different number of atoms by means of Eqs. ( fTTI l and ( fTSl) . As predicted by 
Y. Castin^\ the frequency i2o of the monopole mode (breathing mode) and the 
frequency Q.\ dipole mode (center of mass oscillation) do not depend on A^: 

Q.Q = 2(0 and Qi = (O (23) 

We have found— instead that the frequency Q2 of the quadrupole (/ = 2) mode 
depends on A'^ and on the choice of the gradient coefficient X. In particular, in 
Ref. we have calculated Q2 for increasing values of A'^ up to A'^ = 100. 

In this paper we extend our previous calculations by taking into account much 
larger values of A'^. To excite the quadrupole mode, we solve numerically Eqs. dlTI l 
and ( 1181) with the initial condition 



n(r,f = 0) =«^,(r), (24) 

(25) 



v(r,f = 0) = e(2z^-x''-y^) 
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where n^i(r) is the ground-state density profile and e a small parameter. The re- 
sults are shown in Fig. |3] where we plot the quadrupole frequency D.2 as a function 
of the number A' of atoms for three values of the gradient coefficient X . 

The numerical results of the figure can quite well be captured by a time- 
dependent Gaussian variational approach—. We set 



nix,y,z,t) = 



N 



(26) 



and 



e{x,y,z,t)=l5^{t){x^+y^)+l5,{t)z^ 



(27) 



where G±{t), a,(f), P_L{t), and ji^it) are the time-dependent variational parame- 
ters. Inserting these quantities into the action functional (1211) . after integration over 
space variables we find 



A=N j L{a^{t),a,{t),^^{t)Mt)) dt 
where the effective Lagrangian L is given by 



(28) 



2m 



2o2 



2m 



Of 



2a2 ; 



mO) o 



2^ 



2m 



(aia,)2/3 



(29) 



with g = (3/5)^/2 (37j:2')2/3^ ^ ^ From this effective Lagrangian we calculate ana- 
lytically the small oscillations around the equilibrium configuration, following the 
procedure described in Ref . . In this way we get the following formula for the 
frequency of the quadrupole mode: 



122 = ft) 



2 + 6 



1 + 1 



(30) 



2 gN'-ft 



In the limit A'^ — S- oo it gives the Thomas-Fermi result ^ Q = \f2(i>, while in the 
limit A'^ — > it gives £1 = 2(0, which is the quadrupole oscillation frequency of 
non-interacting atoms ' . The curves reported in Fig. [3] show that this analytical 
formula reproduces remarkably well the numerical results (filled squares). This 
fact suggests the possibility of using Eq. (130b to determine the value of X from 
experimental measurements of 122. 
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6 Collision of resonant Fermi clouds 

We use here the generalized superfluid hydrodynamics equations to obtain the 
long-time dynamics of the collision between two initially separated Fermi clouds, 
simulating the experiments of Ref.^^. The details of our calculations can be found 
m Ref.20. Here we summarize the main results of our study. In our simulations 
we tried to reproduce as closely as possible the experimental conditions of Ref.— 
which we summarize briefly in the following. A 50:50 mixture of the two lowest 
hyperfine states of ^Li (for a total of N = 2x 10^ atoms) is confined by an axially 
symmetric cigar-shaped laser trap, elongated along the z-axis. The resulting trap- 
ping potential is U {r,z) = 0.5m[o)^r^ + 0)?z^], where = +);^, (Or = 2nx 437 
Hz and ft), = 2nx 27.7 Hz. The trapped Fermi cloud is initially bisected by a blue- 
detuned beam which provides a repulsive knife-shaped potential. This potential is 
then suddenly removed, allowing for the two separated parts of the cloud to col- 
lide with each other. The system is then let to evolve for a given hold time t after 
which the trap in the radial direction is removed. The system is allowed to evolve 
for another 1.5 ms during which the gas expands in the r-direction (during this 
extra expansion time, the confining trap frequency along the z-axis is changed to 
ft), = 27r X 20.4 Hz), and finally a (destructive) image of the cloud is taken. The 
process is repeated from the beginning for another different value for the hold time 
f. The main effect observed in the experiment is the presence of shock waves, 
i.e. of regions characterized by large density gradients, in the colliding clouds. 
The experimental results are shown in the right part of Fig. |4] as a sequence of 
one-dimensional profiles obtained by averaging along one transverse direction the 
observed cloud density at different times. 

We simulated the whole procedure by using the Runge-Kutta-Gill fourth-order 
method ^^-^^ to propagate in time the solutions of the following non-linear Schrodinger 
equation (NLSE) 



d 

m—w = 

dt 



+ 2[/(r) +2^(37r2)2/3^ |^'|V3 + (i .4^) ^'l"^! 



4m 2m 4m y¥\ 

(31) 

which is strictly equivalent'^i^ to Eqs. dTTI l and dlSl l, with <f (n, V«) given by Eq. 
(IHl, and 

f(r,f) = V'^^e'^'t'-'') (32) 

Since the confining potential used in the experiments is cigar-shaped, we have 
exploited the resulting cylindrical symmetry of the system by representing the so- 
lution of our NLSE on a 2-dimensional (r, z) grid. During the time evolution of 
our system, when the two clouds start to overlap, many ripples whose wavelength 
is comparable to the interparticle distance are produced in the region of overlap- 
ping densities. In order to properly compare our results with the experimental data 
of resonant fermions , which are characterized by a finite spatial resolution, we 
smooth the calculated profiles at each time / by local averaging the density within 
a space window of ±5 /im centered around the calculated point. 

The results of our simulations, for the whole time duration of the experiments, 
and after the smoothing procedure is applied to the (y-averaged) density profile at 
each time, are shown in the left part of Fig. |4] (see also Ref.^"), plotted along the 
long trap axis, for the same time frames as in the experiment. Remarkably, there 
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Fig. 4 ID density profiles at different times t showing the collision of two strongly interact- 
ing Fermi clouds. Left part: our calculations^". Right part: experimental data from Ref— . The 
normalized density is in units of 10^^/ fim per particle. 



is a striking correspondence between the experimental data and the resuks of our 
simulation. At variance with the current interpretation of the experiments, where 
the role of viscosity is emphasized we obtained a quantitative agreement with 
the experimental observation of the dynamics of the cloud collisions within our 
superfluid effective hydrodynamics approach, where density variations during the 
collision are controlled by a purely dispersive quantum gradient term. 

We find that changing X from the optimal value A = 1 /4 has profound con- 
sequencies on the long-time evolution of the colliding clouds, providing density 
profiles which are completely different from the experimental ones. This is shown 
in Fig|5j where the simulated density profile after 3 ms is calculated using two 
different values for A . Such a strong dependence on A of the time evolution of a 
Fermi cloud made of a large number of atoms is at first sight surprising, because 
the gradient term should become less and less important with increasing A^. We 
believe that such dependence is due to the presence of shock waves (i.e. regions 
characterized by large density gradients) in the colliding clouds, as discussed in 
Ref. 
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Fig. 5 Calculated density profiles after 3 ms. The solid line shows the results obtained using 
X = 0.25, i.e. the same value used to compute the sequence shown in the left part of Fig|4] The 
dotted line shows instead the result obtained using A = 0.1. 

7 Conclusions 

In the first part of the paper we have calculated the density profile and surface 
tension a of the superfluid unitary Fermi gas in a semi-infinite domain by us- 
ing the extended Thomas-Fermi density functional, where the surface effects are 
modelled by the the von Weizsacker gradient term. Indeed we have found that a 
is proportional to VX, where X is the phenomenological coefficient of the von 
Weizsacker term. In the second part of the paper we have investigated dispersive 
dynamical effect which crucially depend on the presence of a gradient term in 
the equation of state of the unitary Fermi gas. In particular, we have studied the 
quadrupole modes of the superfluid system under harmonic confinement, find- 
ing a reliable analytical formula for the oscillation frequency, which reduces to 
the familiar Thomas-Fermi one in the limit of a large number of atoms. Finally, 
we have numerically studied the long-time dynamics of shock waves in the ul- 
tracold unitary Fermi gas. Two main results emerge from our calculations: a) at 
zero temperature the simplest regularization process of the shock is purely dis- 
persive, mediated by the quantum gradient term, which is one of the ingredient in 
our DF approach; b) the quantum gradient term plays an important role not only 
in determining the static density profile of small systems, where surface effects 
are important, but also in the fast dynamics of large systems, where large density 
gradients may arise. 
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